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We calculate second- and fourth-order cumulants of conserved charges in a temperature range 
stretching from the QCD transition region towards the realm of (resummed) perturbation theory. 
We perform lattice simulations with staggered quarks; the continuum extrapolation is based on 
Nt = 10... 24 in the crossover-region and Nt = 8... 16 at higher temperatures. We find that the 
Hadron Resonance Gas model predictions describe the lattice data rather well in the confined phase. 
At high temperatures (above ~250 MeV) we find agreement with the three-loop Hard Thermal Loop 
results. 
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I. INTRODUCTION 

The Quark Gluon Plasma was formed in the Early Universe just a few microseconds after the Big Bang; today it 
is produced in heavy ion collision experiments at the Large Hadron Collider (LHC) at CERN and the Relativistic 
Heavy Ion Collider (RHIC) at Brookhaven Lab (BNL). This phase exists at high temperatures and/or densities, and 
is separated from the hadronic phase of Quantum Chromodynamics (QCD) by a cross-over transition [T]. Lattice 
QCD has determined the temperature of this cross-over in Refs. [2H5]. 

Below the transition temperature, the thermodynamics is governed by massive hadrons with integer charges, whereas 
at high temperature nearly free and nearly massless quarks with fractional charges and gluons dominate. Fluctuations 
of various conserved charges are sensitive probes of the quantum numbers and the associated masses, and have been 
proposed as a signal of the deconfinement transition ISII]. In heavy ion experiments there is an ongoing effort to 
measure the moments of conserved charge distributions [5], which can be related one-to-one to fluctuations. They are 
particularly interesting for the beam energy scan program at RHIC, since they may signal a nearby critical end point: 
higher order moments of net proton distributions are sensitive to an increase in the correlation length . Fluctuations 
can also be used to extract the chemical freeze-out temperature and chemical potential [10] , as an alternative method 
to the thermal fits to particle yields or ratios [IIHIi]- The STAR collaboration has published the first four moments 
of the net-proton m and net-electric charge distributions. In parallel to the experimental efforts, the past years 
have witnessed a rapid development in the lattice calculations of fluctuations [iHl [19] , leading to quantitative estimates 
of the chemical freeze-out temperature and chemical potential for a range of RHIC energies [20] . 

Diagonal quark number susceptibilities have already been calculated in the early dynamical simulations |21H23j : 
these were later complemented by the off-diagonal ones [2^29] . In the following years, higher order fluctuations 
have been calculated up to the sixth order |29L 130] , with the main motivation to extrapolate several thermodynamic 
observables to larger values of the chemical potential. These were staggered simulations projects, but studies with 
Wilson quarks are also emerging [SHSi]. Strangeness fluctuations were used also to locate the transition temperature 
and, for this purpose, they were continuum extrapolated. With Wilson quarks this was done with pion masses down to 
285 MeV [ST] [32] , for staggered quarks the continuum limit was calculated at the physical point [SH!] [36] . Continuum 
results for the other second cumulants appeared first in Ref. m then in Ref. [35]. Selected higher order fluctuations 
were continuum extrapolated first in Refs. [1311351135] . 

Below the cross-over temperature, hadrons (mesons and baryons) dominate the thermodynamics. In this regime, the 
Hadron Resonance Gas (HRG) model provides a simple description of thermodynamic quantities, including specihc 
fluctuations or correlations [niig. Even before simulations with physical quark masses could be performed, lattice 
QCD data were well described by the HRG model if the actual particle spectrum was replaced by the unphysical 
one simulated on the lattice [13| [44] . The success of the HRG model based on the experimental resonance table has 
been demonstrated later in several papers with physical quark masses and continuum extrapolation for the chiral 
condensate [3], the equation of state [l^ and fluctuations [51 [37]. 

The concept of HRG has motivated new studies where fluctuation-based observables were proposed for which, 
within the framework of the HRG model, only particles and resonances with a specific quantum number contribute 
(e.g. baryons in a specihc strangeness sector) [H. Since at low T most lattice results agree with the HRG predictions, 
which is no longer true in the deconhned phase, the highest temperature of agreement can be a model-dependent 
indicator of deconhnement, that can be studied on a havor-specihc basis [33] . 

Very high temperature QCD is best discussed in terms of improved perturbation theory. The QCD thermodynamic 
potential is known up to a^log(a) order [37]. This result was later generalized to hnite chemical potentials [35] and 
the quark number susceptibilities were calculated to the same order |49j . The soft contributions to these unimproved 
perturbative results can be resummed via dimensional reduction |50j . This idea has been applied to the four-loop 
perturbative quark number susceptibilities [SUES]. 

The hard thermal loop perturbation theory reorganizes the perturbative series, enhancing its convergence |53L 154] . 
Recently, the next-to-next-to-leading order pressure and energy density were calculated for the SU(3) theory, [55] . 
dramatically improving the agreement with lattice simulations [Ml [57]. Soon afterwards, the full QCD result was 
calculated, too [55] [53] . Fluctuations were calculated at one [5T] , two m and three loop order |61| . improving the 
earlier HTL calculations of susceptibilities [ST] [53] . This result was later generalized to finite chemical potentials [53] . 

In general, these highly resummed perturbative results are expected to provide a good approximation, but their 
range of validity can only be determined if they are compared with a non-perturbative approach, e.g. lattice QCD 
simulations. Such comparisons have already been made, first on the level of the equation of state [59] . Unfortunately, 
for this observable, the renormalization scale-dependence is too large for a definitive answer on the range of validity. 
Fluctuations, however, offer a more strict test for these diagrammatic approaches because of the rather small renor¬ 
malization scale dependence of the result from dimensional reduction (DR) [52] and from hard thermal loops (HTL) 
m- Today lattice calculations at high temperatures are available e.g. with the HISQ action of the BNL-Bielefeld 
group [351 iO] and also with the 2stout action of the Wuppertal-Budapest collaboration [371 [35]. 
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In this paper, we present results on diagonal and non-diagonal second and fourth order fluctuations, in a temperature 
range which stretches from the transition region to the perturbation theory domain. Our simulations are performed 
within the 2nd generation staggered thermodynamics program (4stout action). We start with the discussion of the 
conserved charges in the grand canonical field theory and provide details on how their fluctuations are calculated on the 
lattice. After describing our lattice thermodynamics program, the scale setting procedure and the finite temperature 
simulations, we highlight the technical challenges of a continuum extrapolation and the estimate of the systematic 
error on the continuum results. The results are organized in two sections. First we consider the cross-over region, 
around the point where the Hadron Resonance Gas loses its predictive power. Afterwards we compare our data to 
(resummed) perturbative results at high temperatures. We close with some concluding remarks pointing to further 
directions of research. 


II. FLUCTUATIONS IN LATTICE QCD 
A. QCD as a grand canonical ensemble 


In a canonical ensemble, the conserved charges are external parameters. In a heavy ion collision, for example, the 
number of baryons, their electric charge and the vanishing strangeness are fixed during the entire collision, expansion 
of the plasma and freeze-out. A grand canonical ensemble emerges if a small sub-system is considered, that is still 
large enough to be close to the thermodynamic limit [65| . 

In QCD there exists a conserved charge for each quark flavor, thus one can introduce four quark chemical potentials 
in a 2 -I- 1 -I- 1 flavor system: /Xs and /ic, in short 

The expectation number of a conserved charge is then found as a derivative with respect to the chemical potential. 

(iv.) = ( 1 ) 

Ufli 


The response of the system to the thermodynamic force fj,i is proportional to the fluctuation of the conserved charge: 




( 2 ) 


Since Ni is an extensive thermodynamic quantity and so is its ^-derivative, there the O(U^) contributions cancel 
in Eq. ([^. Charge conjugation symmetry implies that, at fiq = 0, the expectation value of any odd combination 
vanishes, e.g. the last term in Eq. (i). However, there is no such symmetry for different flavors, allowing e.g. for 
a (NuNd) correlator. The first perturbative diagram that contributes to the latter consists of two fermion loops, 
connected by three gluon lines [62]. 

The free energy density (—T/VlogZ) is proportional to the pressure in large volumes: 

(3) 

The derivatives with respect to the chemical potential can thus be written in terms of the pressure: 


u,d,s,c 


d^+^+k+^{p/T^) 

{dfiuy{dfidy{diis)k{dficy 


(4) 


with fiq = Pq/T. This normalization ensures that the cumulants stay dimensionless, and become finite in the infinite 
volume and infinite temperature limit. In this normalization Xi(T, {/ig}) is the expected number of quarks of the 
given flavor in a volume T~^. 

The higher derivatives with respect to the same quark chemical potential correspond to the higher moments of that 
flavor: 


mean : M = xi variance : = X 2 

skewness: S = X 3 lxl^'^ kurtosis : k = X 4 /X 2 - (5) 

In experiment, the net-charge distribution moments are measured, each carrying an unknown volume factor. A 
known caveat is the fluctuation of these volumes themselves. The study of these goes beyond the scope of this paper, 
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see [SSI EH- For a fixed volume, though, the volume factor can be simply cancelled out by forming ratios of cumulants 
of the same conserved charge: 


ScT = X^/X2 ; no-‘^ = XilX2 

Mla^=X\IX2 ; Sa^ IM = x?,lxi- 


( 6 ) 


Phenomenological models and experiments usually work in the baryon number [B) - electric charge {Q) - strangeness 
(S') basis. Since the charm quark plays a negligible role in the transition region one can express these directions in 
the /X space as a three-dimensional transformation: 


— 

Md = 
^J-s = 


gMS + gMQ : 




AiQ : 




-MQ - MS■ 


(7) 

( 8 ) 

(9) 


The fluctuations of the conserved charges {B, Q and S) can then be expressed in terms of the quark derivatives. 
In addition, the {z component of the) light isospin is often studied with /i/ = (fiu — ^id)■ Assuming zero chemical 
potential and degenerate u and d quarks on the lattice, several simplifications occur, and we have EHEB]: 


X? = ^[2x^+X^+4xn+2xif] , (10) 

X? = ^[5x^+X^-2x)‘r-4xif] , ( 11 ) 

xi = l [x^ - X^il , (12) 

Xi? = l[x^-x^2-xri+x^it] , (13) 

xf/ = -^[X2 + 2Xn] , (14) 

X?f = l[x^2-X^il]- (15) 


Indeed, due to the u d degeneracy the six second order combinations in the B, Q, S space can be expressed in 
terms of four quark correlators. There are 15 fourth order correlators in the (B, Q, S) space that can be expressed in 
terms of 9 fourth order quark-correlators. The kurtosis of the baryon and the electric charge is given by the following 
correlators: 


xf = ^ [2x)( + x! + 6X22 + 12x^1 + 8xil + 8x^r + 8x^i" + 24x^fr + 12x^,fi] , (16) 

X? = ^ [17X4 + Xl + 24x^2^ + 30x^1 - 4x5*1 - 28x^5 - - 24xnl] , (17) 

other, mixed derivatives can be calculated analogously. 

At high temperature, fluctuations approach the Stefan-Boltzmann limit. For an ideal gas, the pressure at finite 
chemical potential reads [551ES] 


p Stt^ 
7^ ^ 



(18) 


For the second and fourth order fluctuations this means that in the high temperature limit X 2 1 and X 4 S/tt^, 
and no mixed derivatives survive. 


B. Fluctuations on the lattice 


The standard way to introduce the chemical potential on the lattice is to modify the temporal links, like the A 4 
component of a homogeneous U(l) field [70]: 


C/4(m) = e^C/ 4 , U+{p)=e-f^U+ 


(19) 
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The fermion matrix M is built from the /r-dependent links. In the staggered formalism, which we will use in this 
paper, each fermion flavor may carry an independent chemical potential. The fermion determinants express a single 
quark flavor. 


Z = J VU e-®»(detM„)i/4(detMd)^/'‘(detM,)i/4(detM,)i/4, 


( 20 ) 


where Sg is the gauge action. To be specific, in this paper we use the tree-level Symanzik improvement in Sg, however 
its form plays no role in the fluctuation-related formulas. The derivative of the staggered fermion matrix M takes the 
following from: 

[Ui{x)^{x + l) + ut{x- b)^{x - 4)] , 

= ]^'no{x) [Ui{x)'ip{x -F 4) - C/ 4 + {x - 6)^p{x - 4)] ; 

any higher odd derivative is equal to dM/d^, while any higher even derivative is equal to d'^M/dfi^. r]y{x) is the 
Kogut-Susskind phase factor. 

For the fourth order /i-derivative one has to evaluate the fourth derivatives of det M. These are traces of the fermion 


matrix that have to be calculated for every generated finite temperature configuration [26] : 

=i|-log(detM,)V4= ItrM-iM', (21) 

Bj = (sgy log(det M,)i/4 = 1 tr (M/M"! - , ( 22 ) 

Cj = ( 3^7 log(det M,) 1/4 = 1 tr (M'M"! - iM] 

-f2MjM-iMjM-iMjM-i) , (23) 

Dj = log(detM,)i/4 = itr {M'^' - iM”M-^M'^' 

+ 12M” M -1 Mj M -1 M'j M- ^ 

, (24) 

Using the simple notation dj for d/dfj.j, the derivatives can now be written for the full free energy: 

d,\ogZ={A,). (25) 


The derivative of the expectation value of any X lattice observable is obtained as 


dg {X) = {XAg) - {X) {Ag) + {Og X) . 


(26) 


When we derive the higher order formulas (see also [26] 1 we assume non-zero chemical potential and use Eq. 
recursively. Setting in the end /i = 0 we have, to second order, 
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d,dg logZ = (AAj) - {A,) {Ag) + 5,g () , (27) 

and to fourth order, exploiting the degeneracy between the light quark flavors: 

df log Z = {A^) - 3 {A^^f + 3 {{Bj) - {B,f) 

+6 {{A^B,) - (Af) {B,}) + 4 (AC,) + (A) , (28) 

d?,dAog Z = (At)-3 (Al)^ 

+3 {(AtB,) - (At) (A)) + (AA) (29) 

dldjlogZ = {At)-3{Al)^ + {Bl)-{BA^ 

+2 {(AtB,) - (At) (A)) (30) 

log Z = (AlA^,) - 2 (A^As)^ - (Al) {A^) + (AA) - (A) (A) 

+ (AlBs) - (Al) (A) + A) - (412) (3^) 







6 


dldAogZ = {AlAP)-?,{Al){A^AP) 

+3((A„Asi?„) — [A^Ag) {BP)) + (AgCu) 

(32) 

dudl log Z = (A„A3) - 3 (A") {A^Ag) 

A-?>{{AuAgBg) — {A^Ag) (Bg)) -I- {A^Cg) 

(33) 

dMl log Z = (A2 A^) - 2 (A„A,)^ - (A^) (A^) 

+ {AlBg)-{Al){Bg) 

(34) 

dldddg log Z = (A3 A,) - 3 (A„A,) {Al) 

+ (A„Asi?„) — (A„As) {Bp ■ 

(35) 


We follow the standard stochastic strategy to calculate the traces A... D, and evaluate them with a large number of 
Gaussian random sources. If one is only interested in up to the fourth derivative, five calls to the linear solver Mx = b 
are necessary for each random source. Since the operator D appears only in connected contributions, we do not need 
it to high accuracy. A, on the other hand, appears in the disconnected term with the most difficult cancellation, so it 
needs to be evaluated more often. A requires one solver, while C requires three solvers. Thus, if we evaluate D with 
N sources, we evaluate the A operator SN times and the B and C operators AN times. 

It was pointed out in [26] that, when products of traces are calculated (e.g. {AA) ~ the two (or more) 

operators in the product must be calculated with different (or uncorrelated) random sources. For this reason, we 
always use quartets of independent sources. We typically use N = 128 quartets in our analysis. Multi-right-hand-side 
solvers are particularly useful in this context, since these typically achieve a higher flop rate on many supercomputers, 
because the gauge fields do not have to be loaded from the memory with each source m- 

The numerical evaluation of these diagrams with multiple random sources can be accelerated by various means. 
One observation was that e.g. the A operator can be split into two parts Aq + 5A, where Aq is the result of a truncated 
solver and 6A is the difference between the truncated result and the full precision solution. The advantage is that 5A 
can be evaluated with less sources, while the more noisy >lo is cheaper to work with m- 


III. LATTICE ACTION AND ENSEMBLES 

This work is part of the second generation thermodynamics program of the Wuppertal-Budapest collaboration. We 
use the tree-level Symanzik gauge action with 2-I-I-1-I flavors of four times stout smeared staggered quarks m, with 
the smearing parameter p = 0.125. 


A. Zero temperature simulations and the line of constant physics 

An essential step, before thermodynamics runs can be started with a new action, is the tuning of the mass parameters 
and the determination of the scale or, in other words, the mass and coupling renormalization of the theory for each 
lattice cut-off that the thermodynamics project intends to use. In this project we use degenerate up and down quarks. 
For simplicity, we do not tune the charm mass separately but accept the continuum extrapolated quark mass ratio 
mc/rris = 11.85 of Ref. m- The light and strange quark masses are obtained by tuning the following ratios to their 
physical values: 


= 27.65, 
J TT 


^phys ^ ^ 1069 

J-rr 


(36) 


where we use the isospin-averaged pion and kaon masses (to^ and tuk ) [IS]. /,r = 130.41 MeV (see Ref. |76|) is used 
to set the scale. 

In this work we use the zero temperature lattice configurations produced for the 4stout T = 0 project m- In the 
lattice spacing range a = 0.188 fm ... 0.077 fm we simulate four or more ensembles for eight inverse bare couplings 
P = 6/g^. The RHMC streams for the ensembles are typically ~ 2000 trajectories long after thermalization. We 
para met rized these ensembles such that they form a ±3% bracket around the physical point, which is defined in 
Eq. (36). The box size of these zero temperature simulations was without exception LtUtt ^ 4. 

In Fig. [^we summarize the zero temperature configurations. For each /3 we interpolated in the space of bare quark 
masses, getting these to a few per mill accuracy. On the left panel of Fig. we show the combinations in Eq. ( |M| . 
The right panel shows the position of individual bare parameters relative to the thus interpolated physical point (with 
details given in Ref. in])- 
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FIG. 1. Summary of our large volume zero temperature runs. The left panel shows the combinations in Eq. ( |36[ l where 
= 1 and Rl/RF^^^ = 1 define the physical point. We determine the physical bare quark masses by interpolating 
the bare parameters to precisely restore the Rs and Rl ratios (see Ref. |77) for details). The right panel shows the bare 
parameters rescaled by the interpolated quark mass. The errors on the rescaled run parameters come solely from the errors of 
the interpolation. 


Our finest large volume ensemble was simulated at /? = 4.0126 on a 96^ x 144 lattice. Its parameters were 
extrapolated and then corrected using simulations at this /3 in the flavor symmetric point, where all three light quark 
masses are degenerate (the charm mass staying physical). 

The tuning effort using the flavor symmetric lattices goes as follows: first, we have to acknowledge that various 
scale setting schemes differ in the cut-off effects. Thus, changing the scale setting or tuning principle may introduce 
different cut-off effects on different parts of the line of constant physics. A continuum extrapolation that spans a 
larger range of lattice spacings will thus be distorted. To prevent this from happening, we match not only the scale 
but also the corrections and check for the insignificance of the effects whenever we are forced to switch between 
scale setting schemes along the line of constant physics. In this particular case, we chose the mass-independent 
renormalization scheme. For a fixed gauge coupling, we define a 3-I-I flavor theory with the bare masses calculated 
from the ones of the 2-I-I-I-I flavor theory: m = |(to„ -I- rud + ms), me = 33.15m. This corresponds to a new scheme, 
and the pseudo-scalar mass to decay constant ratio will have an dependence. We plot this ratio in Fig. (notice 
that, in the 2-I-1+1 theory, m,r// 7 r had no a-dependence by definition). To extract the bare quark masses of the 
2-I-1-I-1 dimensional theory at /3 = 4.00 and /3 = 4.15, we performed several simulations in the 3-1-1 flavor theory and 
interpolated mps//ps in m to match the extrapolation in Fig. We translated the masses back to the 2-I-1-I-1 flavor 
theory. At this point, we had to assume the ms/mu = 27.63 ratio, (which is consistent to our estimate from this 
work) [7^ iTSUSii] . For the large volume simulation at /3 = 4.0126, which was running with such an indirectly tuned 
mass, we show the result in Fig. [l] the physical point is reproduced with an accuracy below one percent. The lattice 
spacings are shown in the plot, for the finest lattice we used the SU(2) low energy constants to extrapolate the final 
one percent to the physical point m- 

For even finer lattices we had to resort to a perturbative continuation of the line of constant physics. For the scale 
setting, the universal two-loop beta function does not yet describe the data. We have an alternative scale setting 
scheme Wq, introduced in (55], which is based on the gradient flow [H3]. In that case, finite volume effects are small 
even for lattices as small as 1.5 fm [52]. This allowed to match again the value and a^-dependence of wq aX /i = 4.1479 
(a « 0.047 fm) and /3 = 4.2562 (a « 0.038 fm). The exploding autocorrelation times have forced us to use extremely 
long update streams (cca. 50000 trajectories) in a 40^ volume. For even finer lattices we again measured and matched 
the flow and its leading lattice artefacts in fixed physical volume and topological sector in several subsequent steps. 
The final scale is plotted in Fig. [^ Since wq is of great interest for a wider community we will discuss its value, 
volume-dependence and other systematics in a publication devoted solely to scale setting. 

Fig. [3] shows two versions of the scale setting. Controlled continuum extrapolations are independent of the choice 
of the scale setting scheme. The equivalence of the schemes on fine lattices is evident from Fig. Nevertheless, this 
choice obviously influences the temperature of a particular ensemble. Especially for observables with large slope in 
temperature (e.g. the quark number susceptibilities in the cross-over region) the scale setting has an impact on the 


















FIG. 2. Using the bare parameters calculated from the 2+1+1 theory’s quark masses m = |(mu + rrid + rris), rria = 33.15m, 
we find a mild dependence for the pseudo-scalar mass-to-decay-constant ratio in the 3+1 flavor (flavor symmetric) theory. 
The matching bare mass m at larger /3 (finer lattice) can be determined at lower computational costs with 3+1 flavors. From 
fti, the bare masses of the 2+1+1 theory can be estimated. 
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FIG. 3. The lattice spacing as a function of the inverse bare gauge coupling. The red squares show the outcome of the 
zero-temperature simulations with Lm,r > 4 for The scale in the wq scheme from the same runs is represented by the red 
circles. The blue dots correspond to smaller volumes, for which we used wo only. The differences coming from the two scale 
setting options are part of our systematic error estimate. 


continuum scaling. We propagate this effect into the final error bars by calculating the continuum limits with both 
scale settings and include this in our study of systematics. 


B. Finite temperature ensembles 

We have generated three sets of ensembles, each with multiple lattice spacings and temperatures. In the first set we 
use the aspect ratio LT = 3, which might have finite volume effects, but gives a more favorable signal/noise ratio than 
larger volumes. The second set has LT = 4 and covers the entire transition range up to 2Tc. Using these ensembles 
we can conclude that, wherever it was possible to perform a meaningful comparison (this includes all second order 
fluctuations and cross-correlators), finite volume effects on the LT = 3 ensembles are negligible for any lattice spacing, 
let alone in the continuum limit which is the largest source of systematic errors. We see significant finite volume effects 
only in the chiral condensate and susceptibility, which are not part of this study. For temperatures T > 300 MeV 
we do not keep the lattice geometry constant in our temperature scan, but keep the physical volume more-or-less 
constant with LTc ^ 2. For the finest, Nt = 16 lattices in this set we have thus used the lattices 80^ x 16, 96^ x 16, 
112^ X 16 and 128^ x 16 for T = 360,440, 520 and 600 MeV, respectively. In the high temperature range, the statistics 
is limited to ~1000 configuration / temperature / lattice spacing. 
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Table [T] shows the statistics for the LT = 4 ensembles in the cross-over region and in the quark gluon plasma phase. 
The temperatures below 150 MeV are used to compare the data to the predictions of the Hadron Resonance Gas 
model. The LT = 3 data set is restricted to the cross-over region (see table 0. In the tables we give the number 
of configurations that we have analyzed for generalized quark number susceptibilities: these are separated by ten 
Rational Hybrid Monte Carlo (RHMC) trajectories. The acceptance range varies between 80 and 95%. 

In the absence of visible finite volume effects in this range, we combine the results of these with the LT = 4 data 
set to enhance the signal. Indeed, the fluctuations of disconnected diagrams (especially — 3(^^) ) in Eq. (28) 
are heavily penalized by large volumes. This contribution also appears in the Taylor coefficients of the expansion 
and is the main source of noise. 


T [MeV] 

32® X 8 

40® X 10 

48® X 12 

64® X 16 

80® X 20 

125 

10514 

10080 

10008 

5027 

2060 

130 

5766 

4625 

10253 

5099 

2000 

135 

14762 

10590 

10060 

10189 

2720 

140 

14863 

5381 

15043 

4959 

5097 

145 

5784 

5020 

10014 

5019 

1280 

150 

5464 

5067 

11043 

5064 

1631 

153 

4985 

5517 

6410 

3641 

- 

155 

5613 

5001 

10137 

5015 

1726 

157 

5526 

5409 

10018 

5160 

1065 

160 

5247 

5017 

4973 

5073 

1082 

165 

8169 

10086 

10496 

5000 

1000 

170 

6005 

6113 

5600 

5111 

600 

175 

12018 

5375 

5058 

5104 

972 

180 

5007 

5089 

5034 

5013 

1000 

190 

4900 

5031 

5121 

5045 

992 

200 

5989 

5002 

6722 

1012 

1000 

220 

5514 

5000 

7231 

1003 

1000 

240 

1712 

5000 

8082 

3947 

1000 

250 

10695 

5685 

5146 

- 

- 

260 

6287 

5000 

8623 

5441 

1000 

270 

11574 

5682 

5684 

- 

- 

280 

7067 

5003 

8751 

1021 

558 

290 

7316 

5680 

5684 

- 

- 

300 

5125 

4917 

5398 

5310 

1011 


TABLE I. The statistics of lattices with LT = 4 aspect ratio. The numbers count the saved and analyzed conhgurations, each 
separated by ten RHMC updates. 


T [MeV] 

24® X 8 

32® X 10 

40® X 12 

48® X 16 

64® X 24 

130 

39161 

7736 

10351 

- 

2007 

135 

41462 

8724 

10696 

9892 

3000 

140 

39867 

8550 

10240 

8248 

1551 

145 

40247 

8518 

10348 

10130 

2550 

150 

39996 

8461 

10569 

6717 

3044 

155 

19953 

8625 

10345 

10211 

1546 

160 

20015 

9174 

11611 

10140 

2063 

165 

10965 

9750 

10219 

10136 

1200 


TABLE 11. The number of analyzed configurations on lattices with LT « 3 aspect ratio. The conhgurations are separated by 
ten RHMC updates. 
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FIG. 4. Fluctuation data from the 4stout action at T = 130 MeV for a lattice resolution range Nt = 10 ... 32. This corresponds 
to the lattice spacings a = 0.15 ... 0.047 fm. A continuum extrapolation is possible by fitting sophisticated models but also via 
a simple linear fit through the finest lattices. The Hadron Resonance Gas model’s prediction (small black arrows) is consistent 
with the continuum limit. 
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FIG. 5. Examples of linear continuum extrapolations for the light-strange correlators Xii and X 2 I at various temperatures. 
The error on our continuum results contains the systematics of varying the fit model, fit window, scale setting and interpolation, 
see text. 


C. Continuum extrapolation 

The continuum extrapolation is mostly based on all available lattice spacings. Since fine lattices have lower statistics, 
the coarsest Nt = 8 results are usually included only in non-linear extrapolations, (e.g. A + B/Nf + CjN^ and other 
variations, where A is the continuum limit). 

While for some observables (e.g. xK^): xiC^)) there is a clear range of safe linear extrapolation (in most cases 
> 10), observables that are related to pion physics (e.g. X?(^); Xii(^)) show a very strong, non-linear 1/Nt 
dependence. Only for very fine lattices (Nt ^ 16) we see a linear regime. Such behaviour have been already reported 
for the second order cumulants IMIET]. 

Here we show the charge fourth and second moment for a single temperature in the confined phase (T = 130 MeV) 
in Fig. 1^ This plot features an additional 96^ x 32 point with 1485 analyzed configurations. We attempt several 
fit models, f\{Nt) = A + B exp{—C/N^) resembles a Boltzmann factor with an artefact mass vanishing as l/N^. 
f 2 {Nt) = A + B/Nt + C/Nt / log{Nt) is similar to including a aa^ term into the extrapolation. The shown continuum 
limit is based on the linear fit only. 

Not all observables require the finest lattices in our data set. Strange quark correlators receive no pion contributions, 
and the small relative taste violation in the kaon sector can be extrapolated away. We find that our data with its 
current precision allow linear htting for Nt > 10. As examples we show the up-strange correlator (xn) and the higher 
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order correlator between the same quarks X 22 Both have only disconnected contributions (see Eqs. ( |27lM| )). 

The parameters of the finite temperature runs have been tuned to have the same temperature in the /„■ scale 
setting scheme. Since we also use the wq scale setting scheme, in that case the temperatures are no longer aligned 
and interpolations are necessary. The alignment of the temperatures is also not perfect in the scale setting scheme, 
thus we interpolate all data sets. The interpolation is performed by fitting a spline through several (7-9) node points 
with two different sets of nodes so that the systematics of the interpolation can be picked up by the systematic error. 
We then perform the continuum extrapolation temperature by temperature, for those temperatures for which we 
had data points. The lattice artefacts of the diagonal fluctuations can be understood from tree-level perturbative 
diagrams |25j . We can correct for the a-independent part of the discretization errors by a T-independent factor (tree 
level improvement) |84] . This factor converges to 1 in the continuum limit. We perform the continuum extrapolation 
in three possible ways: without this improvement, with the tree-level improvement, and with the improvement factor 
of the free energy, that we find empirically to also reduce the cut-off effects at intermediate temperatures. We must 
then judge for every observable separately whether we can include the Nt = 8 and Nt = 10 ensembles, and which 
non-linear models are plausible and match the data. We have given examples for this in Fig. but very often we 
simply add the models A/{1 + B/N^) and A + B/N^ + C/Nf to the linear fit. 


We treat every mentioned option independently and perform 16-32 analyses per temperature, depending on the 
complexity of the continuum scaling. We use this large set of analyses to estimate the systematic errors temperature 
by temperature using the histogram method introduced in Refs. |85L 186) . In this paper we build a histogram of 
the results. The analyses with a fixed data set but different systematics are weighted using the Akaike Information 
Criterion (AIC) [87] . The AIC weigted results corresponding to the various fit windows in 1/N^ are combined with 
uniform weights. In the case of the charm susceptibility we calculate the systematic errors on the finite Nt points 
first and then perform various continuum extrapolations which then enter the histogram method. Since all analyses 
are equally we identify the median with the result. The distribution of results is not necessarily Gaussian and may 
contain isolated combinations of the analysis options that produce outliers. These do not contribute to the median. 
The systematic error is the spread of the distribution. Instead of the standard deviation we use the spread of central 
68 % of the distribution, so that we do not have to make assumptions on the tail of the distribution. The median can 
be calculated for every jackknife or bootstrap sample. We use the variance of the median as statistical error. In the 
plots we show the combined errors, by adding up the systematic and statistical errors in quadrature. 


IV. RESULTS IN THE CROSS-OVER REGION 


Previous works have suggested that the Hadron Resonance Gas (HRG) model provides a good description of the 
data in the range 130-150 MeV jU [351133 Ell SSj , and perhaps missing strange resonances might account for the small 
deviations in the strangeness sector jSSj. 

In this paper we supplement the picture with additional continuum extrapolated data. Finite lattice spacing 
studies (with or without a well improved action) can never state with certainty whether deviations from the model 
are a genuine effect. Here we compare our lattice results using the 2014 edition of the Particle Data Book [89j . 

In our previous paper m we have calculated nearly all the second order fluctuations. Only the most difficult 
correlator was omitted Xii{T), which is not only noisy but had severe lattice spacing effects, similar to xii^) 

Fig. 13 

The continuum extrapolation of Xii{T) and the data in the full lattice spacing range are shown in Fig. together 
with the up-strange correlator Xii{T). The continuum limit for Xif(T) is well described by the HRG model up to 
T « 155 MeV, which lies at the centre of the transition region HIS]. The main hadrons that contribute to the 
HRG prediction are the light mesons, mostly pions (the combination of a quark with an anti-quark makes the Xii 
contribution negative). At high temperatures, heavier hadrons and their resonances have non-negligible Boltzmann 
factors, allowing the baryons (mostly protons) to take over the main role and bend the curve upwards. 

The important role played by the pions is also highlighted by the staggered lattice artefacts (taste breaking) that 
increase the mass of the various staggered pion-like degrees of freedom (tastes) [50] . 

In lattice QGD Xii(F) ~ {A^Ad) /V in Eq. (27)’s notation. The A = (I/4)TrM“^M' operator is a trace over 
the whole lattice. The normalized Gaussian random sources (x) that we use to evaluate A, contribute each as 
X^M“^M'x/4 ^ V. This C-odd estimator is widely oscillating between sources. Thus, in A and then also in the 
stochastic representation of {AA) jV, large cancellations occur between opposite-sign contributions. Refs. [26l]9T] link 
the phase of the fermion determinant at small /is to the odd operators A and C. Indeed, the sign problem is already 
present in the Taylor-expansion technique and in the calculation of baryonic fluctuations in general. 

The consequence is that the severity of the sign problem is related to the magnitude of Xn ■ In early staggered 
studies one saw peak heights of « —0.005 [35], « —0.014 ]35], and « —0.05 are well short of today’s continuum 
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FIG. 6. The up-down correlator (xii) and up-strange correlator (xii) for several lattice spacings and in the continuum limit. 
For the Hadron Resonance Gas model we use the resonance table in the 2014 edition of the Particle Data Book |89| . 


limit in Fig. With the early actions and coarse lattices the calculation of higher derivatives and reweighting were 
easier. 

Note that the light isospin susceptibility (xi) does not depend on the A operator, xi ~ d does not contain 
any disconnected diagrams at all. The fourth derivative xi ~ [6 — {B>)]/V, too, contains only C-even operators. 

Indeed, thermodynamics at finite isospin chemical potential is not plagued by the sign problem. 



100 150 200 250 300 350 400 450 500 


FIG. 7. Light, strange and charm diagonal quark number susceptibilities in the continuum limit as functions of the temperature. 
The quasi-particle model is calculated for a single, non-interacting charm quark with a fitted mass = 1430 MeV. 

A subset of the authors of this paper have remarked that one can observe a hierarchy between flavors in their 
fluctuations |39j . We are now extending the picture and show the continuum extrapolations of the flavor-specihc 
quark number susceptibilities in Fig. The HRG model describes the light flavors reasonably well. The charm 
susceptibility in Fig. ^ rises at higher temperatures, compared to the lighter flavors. It was emphasized in Ref. |92j 
that open charm with fractional baryon charge starts appearing near the chiral crossover temperature. In addition 
to the hadron resonance gas model we show a naive quasiparticle estimate for the charm susceptibility (see also 
|93j). The mass of the charm quark was fitted to the last points = 1430 MeV). This mass is empirical, and 

may depend on the range of the matching to our lattice data. In general the mass of the charm quark is scheme 
dependent. The susceptibility curve runs near the quasiparticle model, qualitatively confirming that X 2 is contributed 
to by the deconfined charm quark. Nevertheless, the quasiparticle model’s results are overestimating the lattice data 
below approx. 350 MeV. This leaves room for multiple interpretations (e.g. T-dependent , limitations of the 
quasiparticle model or charmonium bound states that absorb some of the free quarks). 

Figures]^ and [^detail our continuum results for the fourth order cumulants. The normalized strangeness [39] and 
baryon cumulants HHIISI have been published in earlier works. Here we show the fourth derivative with respect 
to the light single quark chemical potential (Fig. |^. On coarse lattices we see a strong peak around the transition 
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FIG. 8. The single light quark number fourth order susceptibility in the cross-over region. The extrapolation is driven by the 
LT = 3 set of ensembles, which are plotted together with the extrapolation and the HRG prediction. (From T > 165 MeV 
the continuum extrapolation uses only the LT — 4 lattices (not shown)). Since this is a potentially pion-driven observable and 
the Nt = 24 data are not sufficiently precise, the extrapolation is based on A^t = 8... 20 lattices. In the cross-over region we 
consider this a continuum estimate only. 


temperature. 

Such a peak has indeed been expected: if QCD is in the chiral scaling regime with an 0{N) symmetry (i.e. 
the light quark masses are small enough for QCD being nearly chiral) then this scaling is expected to dominate 
the so called magnetic equation of state [94], which parametrizes the singular part of the free energy as a reduced 
temperature and the quark masses that play the role of the magnetic field in the 0{N) model’s language. The chemical 
potential enters through its shifting effect on the transition temperature. At finite fi, the reduced temperature is 
t ^ {T — Tc)/Tc — where k is the curvature of the QCD transition line [95]. Using the critical exponents one 

has, for the n-th derivative, a singular contribution of Xn with 2 — a = PS{1 -I- 1/5) [96]. In the 0(4) 

universality class a = —0.2131(34) [S7]. The non-analytic contribution of is thus singular in the chiral limit 

and has a mild peak near Tc at finite mass, while x^(T) changes sign near Tc |96j . 

The data in Fig.[^show that the peak is strongly reduced on finer lattices, as if we were moving away from the chiral 
limit. It will be interesting to see if this pattern is observed with other actions with an improved dispersion relation. 
Since here the JVt = 24 data have insufficient statistics, we cannot perform a controlled continuum extrapolation at 
all temperatures: we call our result below Tc a continuum estimate. What we see is that already at 145 MeV the 
Hadron Resonance Gas model is unlikely to describe the lattice data. From our extrapolation based on Nt = 8,10,12 
and 16 lattices it is plausible to assume agreement at 135 MeV. 



120 130 140 150 160 170 180 190 200 210 220 


FIG. 9. Continuum limit of the fourth moment of the baryon number distribution. We also show the second moment. 
The HRG model gives the same result for the two observables. The departure of X 4 from X 2 was interpreted as a signal of 
deconfinement in Ref. |46 |. 

The baryon fourth moment shows milder lattice artefacts; here the large statistical errors dominate over the sys- 
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tematic errors (see Fig. |^. We also show x?(^) since the second and fourth moment receive the same prediction 
from the HRG model, independently of how many baryons and mesons are included in the resonance list. The point 
where xf (^) are no longer consistent cannot be described by any resonance list. Multi-baryon states are 

expected to lead to xf > X 2 ^ but here we observe the opposite from T > 155 MeV. The relation yf < X 2 motivates 
the concept that the free energy is dominated by objects with fractional baryon numbers: quarks. Given the trend of 
the HRG model, it is conceivable that the departure point from the HRG model, and the respective maximum of the 
fourth-order derivative {xH or xf) are very close in temperature. 


V. RESULTS AT HIGH TEMPERATURES 


In this section we show our continuum extrapolated results at intermediate and high temperatures. The first 
observables are the off-diagonal quark flavor correlators, already shown in the transition region in Fig. Increasing 
the temperature range (see Fig. 10), we actually see that the value of the light-light correlator spans more than two 
orders of magnitude between Tc and 5Tc. Between ATc and 5Tc the leading perturbative log, which was calculated at 
zero quark mass [62], describes our data. Our data suggest that the light-charm correlator becomes compatible with 
the light-light correlator at about dTj,, but its agreement with the leading log starts a bit earlier. The mass of the 
strange quark is negligible in this observable already at a temperature ~240 MeV. 



FIG. 10. The off-diagonal quark number susceptibilities for various flavor combinations (see also Fig. [^. The light correlator 
spans more than two orders of magnitude in the temperature range between and ST, (using the rescaling factor T, = 
155 MeV). The leading 0{a^ logo) perturbative result is from Ref. [62]. The mass of the strange quark becomes irrelevant 
near 1.5Tc. At STc even the charm quark correlator agrees with the perturbative result, even though the latter was calculated 
at zero mass. 


For the light quark number susceptibility (Fig. 11) there are continuum results available |33I1S|- Here we compare 
to the recent result with the HISQ action (with a combined analysis also using p4 data) |1D|. Our result is compatible 
with both Refs. [ST] |40| within errorbars. Here we also show the latest (improved) perturbative estimates, based on 
hard thermal loops (HTL) jHI] and dimensional reduction (DR) [^. The improvement used in Ref. [S2] has reduced 
the renormalization scale dependence enormously. Our data are approximately one sigma higher than the upper edge 
of the yellow band of the DR result. The central line of the band is calculated at the renormalization scale 27rT, the 
upper edge at 47rT and the lower edge at ttT. 

Both X 4 and Xi are the 


The fourth order cumulants at high and intermediate temperature are shown in Fig. 
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fourth derivative of the free energy with respect to the chemical potential, the difference is that for the former the 
chemical potential is associated with only one of the quarks, whereas for the latter it is associated with all quarks at 
the same time. Here the HTL results have a very small renormalization scale dependence. The data confirms the HTL 
prediction that the Stefan-Boltzmann limit is (almost) reached for xf at intermediate temperatures, x^ approaches 
it much slower. In both cases the improved and resummed perturbative results give an accurate description of lattice 
data above 250 MeV. 

This agreement may seem trivial since the lattice result is continuum extrapolated and resummed perturbation 
theory is evaluated at high temperatures, both approaches are expected to solve QCD. There is a subtle difference, 
however, between HTL theory and lattice solutions. We simulated our ensembles with physical quark masses and 
2-I-1-1-1 dynamical flavors. HTL results, on the other hand, are available for massless quarks only, and for V/ = 3 as 
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FIG. 11. Second order diagonal fluctuations using the single quark chemical potential X 2 vs. the baryon chemical potential 
X 2 i we also compare our data to the BNL-Bielefeld result [40] . 


well as for Nf = 4. The mass of the strange quarks becomes irrelevant before we see agreement between lattice data 
and HTL. At intermediate temperatures the large mass of the charm makes the Nf = 3 Hard Thermal Loop theory 
the closest match to our setting. In order to compare the same observables we do not count the baryon charge of the 
charm quark in xf X 2 ■ To estimate the effect of the charm quark in the sea from the improved perturbation 


theory side we plot the three-flavor and four-flavor result for X 4 together in Fig. 12 (see Ref. |M]b 




FIG. 12. Fourth order cumulants from our lattice study versus hard thermal loops [64] and the result from dimensional 
reduction (DR) ^2]- The small arrows on the right hand side mark the Stefan-Boltzmann limit. 

We close our discussion with the off-diagonal fourth order correlator. In Fig. we show both the light-light and 
the light-strange correlator. Here the effect of the strange quark mass diminishes even sooner, at around 200 MeV. 
The agreement with the HTL result starts at a temperature T ~ 250 MeV, in accordance with the other observables 
We also show the prediction of dimensional reduction [52] . 


VI. CONCLUDING REMARKS 

In this paper we introduced our thermodynamics program with the four-level-smeared (4stout) staggered action. 
We focused on the fluctuations of conserved charges and updated our earlier result on second order fluctuations m- 
Since our first paper on fluctuations, we have introduced very fine lattices (W = 24) in the transition range and 
extended the analysis to high temperatures where a comparison to resummed and improved perturbation theory is 
possible. We have also presented diagonal and off-diagonal fourth order cumulants. Here our data could be used to 
determine the lowest temperature for the three-loop HTL approximation: approx. 250 MeV. 

We have studied whether the hadron resonance gas (HRG) model gives an adequate description of the fluctuation 
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FIG. 13. The off-diagonal fourth order fluctuation at high temperature. Only this off-diagonal derivative has a non-vanishing 
contribution in three-loop HTL [64]. The mass of the strange quark is irrelevant from approx. 200 MeV. Although the 
renormalization scale dependence between ttT and IvrT is large enough to contain the data, an agreement with the central line 
and with its trend in temperature is reached at about 270 MeV. The prediction of the DR method is also shown |52| . there is 
slight disagreement to HTL. Our data is compatible with both at high temperature. 


data. We find that well below the deconfinement temperature, i.e. around 130 MeV, all studied observables are well 
described by the HRG model. This was the most difficult to demonstrate for the fourth moment of the net charge 
distribution Xq, which is a candidate for the freeze-out thermometer at the LHC. In this case, after adding a 96^ x 32 
lattice to the analysis (a = 0.047 fm), our continuum extrapolation based on Nt = 20, 24 and 32 lattices is consistent 
with the HRG model prediction. 

It is very likely that HRG does not describe all aspects of fluctuations in QCD thermodynamics below the transition. 
But for quantities for which it does one can introduce the highest temperature of agreement between lattice and HRG. 
This indicator of deconflnement is unavoidably model-dependent, even if one considers combinations that do not or 
only weakly depend on the actual list of resonances. This temperature can, however, be determined as long as the 
continuum limits are feasible with a sufficient precision. The data on our plots show in most cases an agreement up 
to ~ Tc, which can move to a lower temperature as our precision improves. This should not be confused with the 
limiting temperature of the Hagedorn spectrum, which can be higher. The temperature of highest agreement is not 
the same for all fluctuations as it was also suggested in Ref. [33|, e.g. Xa and very possibly Xa depart from the HRG 
estimates at lower temperatures. This may be a signal of the limitations of the HRG approach, but also suggests that 
the transition is a broad cross-over. 
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